Quantum Monte Carlo calculation of entanglement Renyi entropies for generic quantum systems 
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We present a general scheme for the calculation of the Renyi entropy of a subsystem in quantum many-body 
models that can be efficiently simulated via quantum Monte Carlo. When the simulation is performed at very 
low temperature, the above approach delivers the entanglement Renyi entropy of the subsystem, and it allows to 
explore the crossover to the thermal Renyi entropy as the temperature is increased. We implement this scheme 
explicitly within the Stochastic Series expansion as well as within path-integral Monte Carlo, and apply it to 
quantum spin and quantum rotor models. In the case of quantum spins, we show that relevant models in two 
dimensions with reduced symmetry (XX model or hardcore bosons, transverse-field Ising model at the quantum 
critical point) exhibit an area law for the scaling of the entanglement entropy. 
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Entanglement represents the unique correlation property of 
quantum states, without any classical counterpart, and as such 
it can play a fundamental role in our understanding of quan- 
tum many-body phases from the point of view of non-local 
correlations. The most striking manifestation of entanglement 
in a quantum state \ip) is represented by the mixed nature 
of the reduced density matrix pa describing a subsystem A 
of a quantum many-body system, and defined as the partial 
trace of the total density matrix p — \ij))(xj}\ on the comple- 
ment B, pa = Tr^p. The mixedness of pa can be cap- 
tured by any entropy estimator, the most common being the 
von-Neumann entropy S A = —Tip a logp^, but one can 
equivalently use its generalization, the Renyi entropy (RE) 
CD S^ ] = -log[Tr(p^)]/(l - a), which reduces to von- 
Neumann's in the limit a — > 1. The calculation of entan- 
glement entropies in quantum many-body states appears as a 
formidable task, as it seems to imply the necessity to recon- 
struct the reduced density matrix of a subsystem A; this gen- 
erally represents a hard problem unless A contains very few 
degrees of freedom. In fact this task can be performed effi- 
ciently only in a few cases, including non-interacting bosons 
and fermions on a lattice |2)- For the same models, consid- 
ering fully connected A regions (e.g. hypercubic ones) the 
scaling of the entanglement entropy with the linear size I a 
of the region can be calculated analytically in the asymptotic 
limit Ia — > 00. Analytical and numerical calculations show 



that most models verify a so-called area law S 
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in 

D dimensions, except for critical models in D = 1 and free 
fermions with a (D-l)-dimensional Fermi surface, in which 
the area law is corrected by a multiplicative logarithmic term 
[3 1. Much less is known about models of interacting particles: 
indeed the reduced density matrix can be in principle recon- 
structed efficiently via density-matrix renormalization group 
(DMRG) only in one dimensional lattice systems, while for 
higher-dimensional systems the general method to reconstruct 
Pa is via exact diagonalization, necessarily limited to small 
systems. It is also worth mentioning that the calculation of 
the Renyi entropy for a = 2 has been implemented for SU(2)- 
invariant lattice spin models via a projector Monte Carlo tech- 



nique in Ref . @J Here we propose a new technique to calculate 
the Renyi entropy of a subsystem, valid for arbitrary quantum 
many-body models which admit an efficient quantum Monte 
Carlo (QMC) solution of their equilibrium statistical proper- 
ties. The basic idea is to perform the QMC simulation in an 
extended ensemble for a replicas of the system, treating the 
topology of the (D+l)-dimensional configurations generated 
by QMC as a dynamical variable. We demonstrate this ap- 
proach for the calculation of the a — 2 Renyi entropy both at 
(physically) zero and finite temperature, for two-dimensional 
S = 1/2 quantum spin models with low symmetry, as well as 
for the 0(2) quantum rotor model in D = 1. 

Ref.|5]has shown that, thanks to its trace structure, the Renyi 
entropy at finite temperature can be cast in the form of the log- 
arithm of the ratio of partition functions S a = logR-A /(l — 
a) where R ( « ] = Z A a) /Z a ; here Z a = [Tr (e^ n )] a is 
the ordinary partition function for a replicas of the system, 

(a) 

while Z A is a modified partition function for replicas which 
are "glued" together in the region A. This is best seen in the 
simplest case a = 2, for which 
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Here \riAn B ) is an arbitrary basis of states which are factor- 
ized between the A- and B-region. Regarding e~^ as the 
imaginary-time propagator, Eq. ([T| describes a cyclic propa- 
gation for a time 2f3 of the A-region state, and two indepen- 
dent cyclic propagations for a time f3 of the i?-region state - 
as sketched in Fig. [I] On the other hand Z 2 describes two 
independent propagations for a time 0, 

The statistical mechanics formulation of the Renyi entropy 
has been remarkably exploited in Ref. 5 for the calculation 
of the entanglement entropies for conformal field theories 
(CFT); more recently Refs. 6-8 have implemented a quantum 
Monte Carlo calculation of both Z and Z A separately via 
direct thermodynamic integration of the energy curve E(fi') 
over the interval [0, 0] to obtain the finite-temperature Renyi 
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FIG. 1: Transition from the 2 to the .2^ sector by redefinition 

of the topology of the simulation box in the additional dimension. 

The blue/red regions and lines are associated with the action of the 

imaginary-time propagator e 
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entropy S A (/3). This technique, while being very general, 
appears to be technically limited to finite temperatures, given 
that the statistical error accumulated in the thermodynamic in- 
tegration from infinite temperature down to the temperature of 
interest grows significantly at low T EG). Here we propose 
an alternative QMC approach which cures the above limita- 
tion, allowing to systematically calculate the Renyi entropy of 
a subsystem with a single simulation at the temperature of in- 
terest, performed within an extended ensemble. The central 
idea of our approach is that the ratio of two partition functions 
can be generally estimated with Monte Carlo by performing 
a simulation in an ensemble which is the union of the two, 
Z 2 U Z A . Whichever quantum Monte Carlo approach is 
used for the estimation of the equilibrium statistical proper- 
ties of the Hamiltonian %, it should allow to write Z A ' in the 
form 
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where C = (tia, Ub', m-A, m B\1-') is a QMC configuration, 
in which the state \ua) is propagated to \itia) and then back 
to itself in the region A, while the states \ub) and |m.e) are 
propagated onto themselves independently, and the propaga- 
tion scheme is represented by V: V is generally a path in the 
computational basis \i/j(t)) = \iI>a(t))\i(>b(t)) parametrized 
by the (continuous) imaginary time r £ [0, 0\ as in path- 
integral Monte Carlo (PIMC) [9], or by the propagation step 
index r = p, associated with a string of bond operators, as 
in Stochastic Series Expansion (SSE) iflOl . Within the above 
notation, Z 2 = Z A ^ . wa is the statistical weight of a con- 
figuration; the Hamiltonian H lends itself to an efficient QMC 



simulation if wa > for all configurations, and if the weights 

wa can be calculated efficiently. 

Our method is then based on constructing a simulation 

which moves dynamically between the Z 2 ensemble and 

(2) 
the Z A while respecting detailed balance condition. The 

move from one ensemble to another can be performed with 

Metropolis probability 



P IZ' S -> Z\ 
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w A {C) 
wa=sz{C) 
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and viceversa for the reverse move. The partition function 
ratio R A ' is then simply estimated as 



p (2) _ / N A 
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where Na is the number of MC steps in the ensemble with a 
given region A, and (...)mc is the Monte Carlo average. A 
straightforward generalization of the above formulas is pos- 
sible for a > 2. We would like to point out that an analo- 
gous extended-ensemble QMC scheme is the one defining the 
QMC estimator for observables which are off-diagonal in the 
computational basis [11]. 

In practice, for a simulation on discrete degrees of free- 
dom on a lattice - e.g. quantum spins, lattice gases - the 
weights wa(C) and Wa=0 {C) cannot be simultaneously non- 
vanishing unless the condition \ua) = [tUa) is satisfied, in 
which case the transition in the propagation topology (Fig.[T| 
is microcanonical, namely wa(C) = wa=z{C). In the case 
of continuous lattice variables - e.g. quantum rotors - or of 
particles in continuous space, the "rewiring" of worldlines 
demanded by the transition between the ensembles can be 
in principle always performed, although it will have an ac- 
ceptance rate which is low if the configurations \ha) and 
\itia) are very different. Assuming to use a PIMC scheme 
in which the imaginary time is discretized in steps At, and 
indicating with "Hod the part of the Hamiltonian which is 
off-diagonal in the computational basis, one has that the ra- 
tio wa(C)/wa=0(C) takes the expression 
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As indicated in Fig. [T] the primed configurations are those 
which are connected to the un-primed ones by a single propa- 
gation step. 

The above scheme provides an efficient estimate of R A , 
and therefore of S A , by performing a single simulation at 
the temperature of interest. When the temperature is chosen 
to be so low as to remove thermal effects on a finite-size sim- 
ulation box, one can gain access to the entanglement Renyi 
entropy. Unless otherwise specified, in the following we will 
show simulation results for the general case of an XYZ Hamil- 
tonian in a field 
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where Sf are S = 1/2 spin operators, J > 0, and (ij) in- 
dicates a pair of nearest neighbors on a D-dimensional hy- 
percubic lattice. A validation of our approach comes from 
the comparison with exact results in D = 1, which are avail- 
able e.g. for the case of the XX model {A y = 1, A z = 0); 
such a model admits a mapping onto a system of free fermions 
|[T2l . whose entanglement properties can be calculated from 
the knowledge of two-point correlations [2|. Fig. [2ja) shows 
the data for a L = 64 chain with periodic boundary condi- 
tions, simulated with the SSE algorithm; very good agree- 
ment is found between the exact results and QMC results at 
an inverse temperature [3 J = 200. In principle the whole 
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FIG. 2: Upper panel: Entanglement entropy of the ID XX chain; 
the solid line corresponds to exact diagonalization. Lower panel: 
Entanglement entropy of the ID 0(2) quantum rotor model. Here 
L — 64, t = k B T/J, and e - Aril; increments AZ - 1 
are used; the dashed lines correspond to fits to the CFT prediction 
(l/4)\o g [C(l A \L)]+ Sl . 
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formally rewriting R A as 
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where Ai is a sequence of N blocks of increasing size such 



that A Q 



and A 



N 



A. Each of the ratios R 
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A it Ai 



can be 



estimated efficiently, as it represents the ratio between the par- 
tition functions of systems which are 2/3-periodic on regions 
Ai and Ai+i chosen so as to differ only by a few sites (or by 
a few interparticle spacings in continuum space). The Renyi 
entropy is then the sum of contributions from the successive 



increments AAi that lead from to A, S A — Yl S 



where S 
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ments AAi gi ve rise to inefficient estimates of the correspond- 
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too small ones lead to a sizable accumu- 



ing ratios R AijA(+1 , 
lated error on the sum; yet an optimal size of the increment 
can be found minimizing the final error on S A . In Fig. 12(a) 
we have used linear increments of size AAi = Al = 5. 
Nonetheless it can still be seen that the precision of the results 
is not optimal for I a > 20. This is a result of the slow in- 
crease of entanglement entropy in ID systems (especially for 
I A ~ L/2): if the ratios i?( 2 ' are known with a given relative 
error cr = AR^ / R^ 2 \ the corresponding entanglement in- 
crement S^ 2 ) has a relative error es = £r/S^ 2 \ which can be 
much bigger than en when the increment is small. This means 
that the QMC technique enjoys a faster scaling of the entropy, 
as found e.g. in 2D systems or at finite temperature; as we 
will see, the quality of the 2D data is significantly better. 
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S A curve as a function of I a can be obtained by perform- 
ing simulations in the joint Z 2 U Z A ensemble; in prac- 
tice, nonetheless, the transition rate between the two ensem- 
bles is strongly suppressed when the size of A grows, given 

that the condition \n A ) = I'm a) is increasingly hard to sat- 

(2) 
isfy. This aspect reflects the fact that the ratio R A esti- 

2(2) 



mated in the simulation decreases exponentially with S A , 
In the D = 1 case in question, in 
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the case in which an area law holds, one has an even more se- 



rious decrease R 
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exp(- 



-b I a l )- Hence the events that 



the simulation has to count become increasingly rare, so that 
the simulation length should naively scale as (R 
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To 



cure this problem we use the increment trick from Ref. [4] by 
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FIG. 3: Entanglement entropy 2D spin models with various sym- 
metries; the error bars are smaller than the symbol sizes. The dashed 
lines represent fits to the equation in the main text. 

Having validated the approach against exact results, we can 
apply it to yet unexplored models. To demonstrate the versa- 
tility of the QMC RE estimator, we apply it to the study of 
a model with continuous quantum variables, namely the ID 
0(2) quantum rotor model H = — 2J^/- , cos((j>i — cf>j) — 

^ ^ li d 2 /d<p 2 , in which <fn G [0, 2tt\. This model represents 
an approximation to the Bose-Hubbard model with hopping 
J and repulsion U for large integer filling, and it exhibits a 
superfluid-insulator quantum phase transition for increasing 



U/J. Such a model can be studied via PIMC [ 1 3 1 with dis- 
cretized imaginary time (in steps At). Fig. |2|b) shows the 
RE for a chain of L = 64 sites at variable U/J; we ob- 
serve that for sufficiently small U/J the RE obeys the CFT 
prediction S A 2) = (c/4)log[C(l A \L)] + s 1 where c = 1, 
C(x\L) = L/TTsm(Trx/L), and S\ = s\{U / J) is a constant 
dependent on the Hamiltonian parameters. On the other hand 
for large U/J the CFT prediction is no longer verified, as the 
system enters an insulating gapped phase. 

We then move to 2D systems, and consider three repre- 
sentatives of the three symmetry sectors of the XYZ model in 
zero field, namely the case of a SU(2) invariant Heisenberg (or 
XXX) model (A y = A z = 1), the case of a U(l) symmetric 
XX model {A y = 1, A z = 0), and the case of the Z2 symmet- 
ric anisotropic XY model (A y = 0.8, A z = 0). In all three 
cases we consider A regions with a square geometry I a x I a, 
grown in linear increments of (typically) 5 sites, and we plot 
the data as a function of the region boundary I = 4(Ia — 1). 
The simulations have been performed with the SSE algorithm 
on lattices with L x L size up to L = 36, and at a temperature 
j3J m L guaranteeing the removal of thermal contributions. 
The case of the XXX model has been previously investigated 
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FIG. 4: Renyi entropy for the 2D XX model at increasing tempera- 
ture. Error bars are smaller than the symbol size. 

in Refs. |01[T4] via projector QMC, and we confirm their find- 
ing of an area law scaling of entanglement entropy. As clearly 
shown in Fig. [3] an area law is also observed for the other two 
models with reduced symmetry: in all three cases the scaling 
of the Renyi entropy is very well fitted by an area law plus 
subleading corrections, /(/) = b I + clog I + d; the fit coeffi- 
cients (obtained by discarding data with I < l m [ n ) are reported 
in Table II] [T5]. In particular, the coefficient b of the domi- 
nant area-law term decreases systematically as the symmetry 
of the model is decreased; this is consistent with the picture 
that a lower symmetry confines quantum fluctuations to a re- 
stricted region of spin space, thereby lowering entanglement 
properties. 

The 2D XX model maps onto hardcore bosons, and it is 
directly relevant e.g. to current cold-atom experiments. To 
make contact with a more realistic experimental situation, we 
have studied the effect of an increasing temperature on the 
scaling of the Renyi entropy, an aspect which can be quite 



model 


b 


c 


d 


''min 


2D XXX (A y = A z = 1) 


0.099(1) 


0.48(4) 


0.05(10) 


16 


2D XX (A s = 1, A z = 0) 


0.045(1) 


0.31(3) 


0.52(5) 


8 


2D XY (A y = 0.8, A z = 0) 


0.013(1) 


-0.02(2) 


0.76(2) 


8 


2D TFI (A y = A z = 0) - QCP 


0.0332(4) 


-0.03(1) 


0.15(2) 


8 



TABLE I: Fit coefficients for the three models investigated in Fig. [3] 



naturally investigated with finite-T QMC. Given that thermal 
entropies are extensive, on general grounds one expects finite 
temperatures to introduce a volume law in the scaling, namely 
a a(T)l A term fatally masking the area-law term. Nonetheless 
the growth of the volume-law coefficient with temperature, 
a(T), appears fairly slow: as shown in Fig. Han area law is 
also observed at a moderate, finite temperature, T/J = 0.1, 
for the block sizes considered here (the biggest being I a = 
16), implying that a(T)l% < Ab(l A - 1). On the other hand, 
at a temperature T/J = 0.5 and higher the volume law term 
dominates already for small block sizes. These results point 
at the fact that area laws of Renyi entropy are observable even 
at finite temperature and for moderate block sizes, which are 
indeed relevant for experiments. 

We conclude our discussion of D — 2 models with the 
case of the 2D transverse-field Ising (TFI) model, A y = 
A z = 0, which displays a quantum critical point (QCP) at 
H c sa 1.52J |[T6l . Our approach enables us to investigate 
this two-dimensional quantum critical system in search for 
special entanglement signatures. As shown in Fig. [3] and in 
Table II] we observe that the entanglement RE obeys an area 
law with a negative logarithmic correction and a positive ad- 
ditive constant. These findings are in quantitative agreement 
with recent field-theory results for a QCP with dynamical crit- 
ical exponent z = 1, predicting universal negative logarithmic 
corrections coming from corners of the A region |[T5l[T7l . and 
a universal (a-dependent) additive constant in D = 2 ifTHl . 
This shows that QMC simulations can quantitatively extract 
the subleading corrections; their universality can be directly 
tested by investigating different microscopic models exhibit- 
ing QCPs in the same universality class. 

In conclusion we have demonstrated a simple approach to 
incorporate an estimator of subsystem Renyi entropies into 
any finite-temperature quantum Monte Carlo scheme. This 
approach complements the estimator developed within the 
projector-QMC scheme ||4), and it paves the way for a sys- 
tematic investigation of entanglement entropies in a large va- 
riety of interacting quantum systems in arbitrary dimensions, 
such as quantum fluids, quantum spin systems, O(N) quan- 
tum rotor models, quantum field theories etc., as long as they 
are accessible to a QMC study. Contrary to DMRG or to vari- 
ational methods based on tensor-network states, QMC sim- 
ulations are completely unbiased with respect to the scaling 
of entanglement, and the QMC entanglement estimator actu- 
ally performs better the faster the entanglement grows with the 
subsystem size. Moreover QMC represents a natural platform 
to investigate the statistics of local quantum fluctuations in re- 



alistic systems, in the attempt to relate measurable fluctuation 
properties with entanglement properties [ 19 1. 

T. R. acknowledges fruitful discussions with J. I. Cirac and 
R. Melko which have sparked the present study. 
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SUPPLEMENTARY MATERIAL 

Here we provide a more detailed discussion of the co- 
efficients of the Renyi entropy scaling of a square subsys- 
tem embedded in a 2D quantum spin system. We fit our 
QMC data for 2D quantum spin models to the form f(l) = 
bl + c log I + d. On a L x L lattice, the fits are performed over 
a region of boundary sizes [7 m j n ,4(L/2 — 1)] whose lower 
bound l m i U is gradually grown to check convergence. Fig. [5] 
shows the fit coefficients as a function of Z min . As a general 
criterion, if convergence of the fitting parameters within the 
error bar is achieved for a given ^ in , we choose as best fit pa- 
rameters (shown in Tablell} the ones corresponding to l* nin — 4 
(given that they are consistent with the Z* nin values and have 
smaller error bars). In general we observe that the coefficient 
of the area law b is very stable to variations of Z m i n , while 
shrinking too much the fitting region leads to a transition in 



the coefficients of the subleading terms c and d. Nonethe- 
less we observe that convergence in the fitting coefficients is 
achieved before the transition; we observe that the transitions 
are systematically accompanied by a degradation in the pre- 
cision of the resulting fit coefficients, and we argue that they 
can be attributed to the limited data sets that are left to fit if 
Zmin grows too big - such limited data sets can only provide 
reliably the coefficient of the leading term, but hardly those of 
the subleading ones. 
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FIG. 5: Evolution of the fit coefficients of S\ when increasing the 
lower bound of the fit region [l m i n , 4(L/2 — 1)]. 

Recent field-theoretical studies DHU have pointed out that 
several two-dimensional quantum systems should exhibit a 
dominant area-law scaling of the entanglement entropy with a 
non-universal coefficient b - given that such coefficient would 
depend on the short-distance cutoff related to the details of 
the microscopic model of origin. On the other hand, the sub- 
leading logarithmic and constant terms can indeed be cutoff 
independent and universal. In the case of z = 1 quantum crit- 
ical points (QCPs) 1 3 1, if the A region has a smooth boundary 
one expects that 6 = and an additive universal constant d, 
dependent on the index a of the considered Renyi entropy. 
Our finding for the QCP of the 2D transverse field Ising (TFI) 
model is that d is indeed finite and positive for a — 2; in 
the case a = 1 (von-Neumann's entropy) Ref. 5 estimates as 
well a positive additive constant. On the other hand, we also 
find a finite logarithmic correction with a negative coefficient, 
c = —0.03(1). This is not at all surprising, given that the long 
wavelength properties of the QCP in the 2D TFI model should 
be represented by a free relativistic theory. For such a theory 
Ref. [H predicts that a region A with corners on the boundary 
will acquire universal negative logarithmic contributions; in 
particular each comer should contribute a term s=s —0.0062 to 
the c coefficient, which for the 4 corners of square A regions 
provides a value in quantitative agreement with our estimate 
of the c coefficient. 

If the ground state has a finite correlation length (as it is 



the case for the 2D anisotropic XY model), one expects a 
correlation-length dependent additive constant d |3); the ex- 
isting predictions of logarithmic corrections do not apply to 
this case. Indeed we find a sizable additive constant, and a 
logarithmic correction which is consistent with zero. Finally, 
in the case of the 2D XX and XXX model, the ground state 
has an infinite correlation length and it develops long-range 
order in the thermodynamic limit. For the case of the 2D 
XXX model, fits to projector QMC data have been performed 
in Ref. 6 using the function /'(7a) = 4a'l A + c' log(4^) + d' 
of the boundary size estimated as 41a (this estimate double- 
counts the corner spins). When fitted to the /' function our 
data deliver coefficients which are indeed in agreement with 
the ones quoted in Ref.|6l 

Both the XX and the XXX model have linearly dispers- 
ing gapless Goldstone modes (two for the XXX model, and 
one for the XX model), each described in the long-wavelength 
limit by a free relativistic theory. Following Ref.|T|one would 
expect negative logarithmic corrections coming from corners, 
but in fact our results point at a positive c coefficient for both 
models. This result is consistent with what was initially found 
numerically in Ref. 6 where positive logarithmic corrections 
have been shown to exist for the XXX model even in absence 
of corners. Prompted by the results of Ref. |6] Ref. @]has re- 
cently predicted that in systems exhibiting spontaneous sym- 
metry breaking in the thermodynamic limit, one should expect 
a positive logarithmic correction with a coefficient c which 
takes the simple form Nq(D — l)/2 where Nq is the num- 
ber of Goldstone modes. This would imply that c = 1 for 
the 2D XXX model and c = 1/2 for the 2D XX model. Our 
observation is not consistent with this prediction, even when 



taking into account negative logarithmic contributions com- 
ing from corners. Nonetheless we observe that the coefficient 
of the XXX model is significantly larger than that for the XX 
model. One might argue that a possible source of discrep- 
ancy between our results and those of Ref. [4] stems from the 
finite-temperature nature of our data. Indeed a temperature 
T ~ L^ 1 (at which our simulations are conducted) is suffi- 
ciently low to eliminate Goldstone-mode excitations, but not 
to eliminate the thermal occupation of the low-lying tower- 
of-states excitations [7 1 (to eliminate those states one would 
need a prohibitively low temperature, T ~ L~ 2 ); yet Ref. 4 
suggests that in the case L~ 2 <C T -^i L^ 1 their prediction 
should still hold. A further source of discrepancy could be 
the fact that our finite-size results fail to correctly capture the 
behavior of the subleading terms in the limit / — > oo. Future 
larger-scale simulations should be able to clarify this issue. 
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